Load data

Intensity and shapefile data is handled here.

Make Thresholded Raster for 2017, Same Threshold as 2018

Cut 2017 VV TIF to the used extent.

allRas2017 <- raster("allRas2017.tif")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
allRas2018 <- raster("allRas2018.tif")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
colo = viridisLite::inferno(29)
breaks = seq(0, 29, 1)
image(allRas2017, col = colo, breaks = breaks, main="2017")
colo = viridisLite::inferno(31)
breaks = seq(0, 31, 1)
image(allRas2018, col = colo, breaks = breaks, main="2018")

change <- raster("norm_change_17_18.tif")
colo = viridisLite::inferno(20)
breaks = seq(-1, 1, 0.1)
image(change, col = colo, breaks = breaks, main="Normalized Change 2017 - 2018", xlim=c(730000, 770000))

plot2017 <- loadPolygon(shape[124,], dates2017)
plot(plot2017[1,,,])

Rain 2017

## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Lets See What Happened End of June

# plot overview
area1 <- loadPolygon(study_area[1,], dates2017)
plot(area1[,,,12:17])

The second big precipitation event is 24/07/2019

# plot overview
# area1 <- loadPolygon(study_area[1,], dates2017)
plot(area1[,,,16:17])

before <- as(area1[1,,,16], "Raster")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): Discarded
## datum European_Terrestrial_Reference_System_1989 in CRS definition
after <- as(area1[1,,,17], "Raster")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum European_Terrestrial_Reference_System_1989 in CRS definition
# see histograms to evaluate if the scene has gotten brighter
hist(before)

hist(after)

Plot Scenes with Threshold, First Prec Event

threshold <- -15.34187
image(area1[,,,14], col = c("white", "black"), breaks = c(-35, threshold, 35), main="25.06., White is Smaller -15.3")
image(area1[,,,15], col = c("white", "black"), breaks = c(-35, threshold, 35), main="07.07., White is Smaller -15.3")

threshold <- -20.49962
image(area1[,,,14], col = c("white", "black"), breaks = c(-35, threshold, 35), main="25.06., White is Smaller -20.5")
image(area1[,,,15], col = c("white", "black"), breaks = c(-35, threshold, 35), main="07.07., White is Smaller -20.5")

A gain in values below -20 is observed, together with a reduction of surfaces below -15. This could be a hint that -15 is not a usable water threshold.

See as RGB (even better in QGis)

stack <- raster::stack(as(area1[1,,,13], "Raster"), as(area1[1,,,14], "Raster"), as(area1[1,,,15], "Raster"))
plotRGB(stack, r=1, g=2, b=3, stretch="lin")

Bigger Area, Bring Back the LazyLoad Function

Which however only works when giving more than 1 time step which was useless when I wrote it but is perfect here. It extracts the requested timesteps fom the stars_proxy objects and therefore reduces computing time and resources significantly, allowing for larger areas to be loaded.

area2 <- lazyLoadPolygon(study_area[2,], c(13:15))
stack <- raster::stack(as(area2[1,,,1], "Raster"), as(area2[1,,,2], "Raster"), as(area2[1,,,3], "Raster"))
plotRGB(stack, r=1, g=2, b=3, stretch="lin", main="RGB of 13.06., 25.06. 07.07.2017")

New Water Threshold Training with Reworked Water Shape

# prepare data
val.svm$type <- as.factor(val.svm$type)
attach(val.svm)
x <- subset(val.svm, select=-type)
y <- type
# make model
svmmod <- svm(x,y)
# summary(svmmod)
# make prediction, confusion matrix
pred <- predict(svmmod,x)
table(pred, y)
##        y
## pred    field water
##   field   773    12
##   water     0   264
# threshold is
svmmod$x.scale$`scaled:center`
## [1] -16.38373
# plot classified with threshold calculated by SVM
plot(area1[1,,,1], col = c("white", "black"), breaks = c(-35, svmmod$x.scale$`scaled:center`, 10))
# vs old threshold
# plot(all[1,,,1], col = c("white", "black"), breaks = c(-35, thresh$threshold, 10))
ggplot(val.svm, aes(x=type, y=VV)) +
  geom_boxplot() +
  geom_hline(yintercept = svmmod$x.scale$`scaled:center`)

Look at the End of June Again

plot(area1[,,,14])
plot(area1[,,,15])
threshold <- -15.34187
image(area1[,,,14], col = c("white", "black"), breaks = c(-35, threshold, 35), main="25.06., White is Smaller -15.3")
image(area1[,,,15], col = c("white", "black"), breaks = c(-35, threshold, 35), main="07.07., White is Smaller -15.3")
threshold <- svmmod$x.scale$`scaled:center`
image(area1[,,,14], col = c("white", "black"), breaks = c(-35, threshold, 35), main="25.06., White is Smaller -18")
image(area1[,,,15], col = c("white", "black"), breaks = c(-35, threshold, 35), main="07.07., White is Smaller -18")

Next Steps

Sneak Peak: 2017 2018